 %%%%%%%%%%%%%%%%QPSK%%%%%%%%%%%%%%%%

clear; clc;
%liczba bitów stosowanych w tecie
nbit=20;
%wartociowoæ modulacji
M=4;
L=sqrt(M);

%How many pathways (Tx x Rx)
N=4;


%Entityattributevalue model
Eav=4;
%Input signal from 0 to 3
x_signal_input = randint(nbit,1,4);
%k parameter needed for normalization
k=sqrt((3*Eav)/(2*(M-1)));

%macierz wartoci nadawanych
Tx=[];

%Making of QPSK constelation
for i=1:nbit
    if x_signal_input(i,1) == 0
        Tx=[Tx (k+j*k)];
    elseif x_signal_input(i,1) == 1
        Tx=[Tx (-k+j*k)];
    elseif x_signal_input(i,1) == 2
        Tx=[Tx (k-j*k)];
    else
        Tx=[Tx (-k-j*k)];
    end;
end;
numel(Tx)

figure();
plot(real(Tx),imag(Tx),'ko','MarkerFaceColor','b');
title('QPSK INPUT');

%Gaussian noise
% 0 - noise power
% 1 - max value
y2 = wgn(nbit,1,0,'complex');

%RAND NOSIE
noise=randn(1,4*nbit)+j*randn(1,4*nbit); % 4* bo kazdy symbol jest przesylany 2 razy przez 2 drogi
%SIGNAL TO NOISE RATIO
SNR  = 10;

No = 10^(-SNR/10);

%Adding scaled noise to input signal
%Tx = repmat(Tx,1,2); % kazdy 
%Rx=Tx+sqrt(No/2).*noise;

%Dodanie do wektora nadanych symboli powtorzonych w chwili t1
%w momencie t1 antena 1. nadaje odwrocone sprzezenie sygnalu nadanego przez
%antene 2. w momencie t0
Tx_second_antenna = zeros(1,length(Tx));
for i=1:2:(length(Tx)-1)
    Tx_second_antenna(i)=-Tx(i+1)';
end
%w momencie t1 antena 2. nadaje sprzezenie sygnalu nadanego przez antene 1.
%w momencie t0
for i=2:2:(length(Tx))
    Tx_second_antenna(i)=Tx(i-1)';
end
%Tx=[Tx, Tx_second_antenna];
%kazdy symbol przebywa 2 drogi - *2
%Tx=repmat(Tx,1,2);

%dodanie szumu
%Rx=Tx+sqrt(No/2).*noise;

%Apply Rice channel propagation model
%sigma = std(Tx) ; %scattered mean power
%A = mean(Tx); %Amplituda składowej dominującej
%A=Tx;
%I0=besseli(0,(A.*Rx/(sigma).^2));
%Rx=(Rx/(sigma).^2) .* exp(-((Rx).^2 + A.^2)/(2*sigma.^2)) .* I0;
%Rx=(Rx/(sigma).^2) .* exp((-(Rx).^2 + A.^2)./(2*sigma.^2)) .* I0;

%sygnaly odbierane - zbedne raczej
h11=zeros(1,nbit);
h12=zeros(1,nbit);
h21=zeros(1,nbit);
h22=zeros(1,nbit);
y1=zeros(1,nbit);
y2=zeros(1,nbit);
y1p=zeros(1,nbit);
y2p=zeros(1,nbit);
clear i;
for j = 1:nbit
    h11(j)=ricernd(0.5,0.25)+i*ricernd(0.5,0.25);
    h12(j)=ricernd(0.5,0.25)+i*ricernd(0.5,0.25);
    h21(j)=ricernd(0.5,0.25)+i*ricernd(0.5,0.25);
    h22(j)=ricernd(0.5,0.25)+i*ricernd(0.5,0.25);
    %y1(j)= h11*Tx(i)+h12*Tx(i+1) + noise(i);
    %y2(i)= h21*Tx(i)+h22*Tx(i+1) + noise(nbit+i);
    %y1p=h11*(-Tx(i+1)') + h21*(Tx(i)') + noise(2*nbit + i);
    %y2p=h21*(-Tx(i+1)') + h22*(Tx(i)') + noise(3*nbit +i);
end

xodeb=zeros(1,nbit);
%odbior Alamoutiego - tu trzeba gdzies chyba h11 itd z modelu rice'a wziac?
for i=1:2:nbit-1
    xodeb(i)=(Tx(i)*((abs(h11(i)))^2 + (abs(h12(i)))^2 + (abs(h21(i)))^2 + (abs(h22(i)))^2) + h11(i)'*noise(i) + h12(i)*noise(nbit+i)' + h12(i)'*noise(2*nbit+i) + h22(i)*noise(3*nbit+i)')/((abs(h11(i)))^2 + (abs(h12(i)))^2 + (abs(h21(i)))^2 + (abs(h22(i)))^2); %38 strona wyklad
end
for i=2:2:nbit
    xodeb(i)=(Tx(i)*((abs(h11(i)))^2 + (abs(h12(i)))^2 + (abs(h21(i)))^2 + (abs(h22(i)))^2) + h11(i)*noise(i)' + h12(i)'*noise(nbit+i) + h12(i)*noise(2*nbit+i)' + h22(i)'*noise(3*nbit+i))/((abs(h11(i)))^2 + (abs(h12(i)))^2 + (abs(h21(i)))^2 + (abs(h22(i)))^2); %38 strona wyklad
end
Rx=xodeb;


%OUTPUT RES
figure();
plot(real(Rx),imag(Rx),'ko','MarkerFaceColor','b');
title('QPSK OUTPUT');



sym1=zeros(1,length(Rx));
sym2=zeros(1,length(Rx));
sym3=zeros(1,length(Rx));
sym4=zeros(1,length(Rx));


for i=1:(length(Rx)/N)
    for j=0:(N-1)
        if(real(Rx(i+j*(length(Rx)/N)))<0 ) %&& imag(Rx(i+j*(length(Rx)/4))>0))
            if(imag(Rx(i+j*(length(Rx)/N)))>0)
                sym1(1,i)=sym1(1,i)+1;
            else
                sym3(1,i)=sym3(1,i)+1;
            end
        end
        if(real(Rx(i+j*(length(Rx)/N)))>0 )%&& imag(Rx(i+j*(length(Rx)/4))>1))
            if(imag(Rx(i+j*(length(Rx)/N)))>0)
                sym2(1,i)=sym2(1,i)+1;
            else
                sym4(1,i)=sym4(1,i)+1;
            end
        end
    end
end

%M=length(Rx)/N;
M=nbit;
recsym=zeros(1,nbit);

numel(Rx)
numel(sym1)
for i = 1 : nbit
    if(sym1(i)>=sym2(i) && sym1(i)>=sym3(i) && sym1(i)>=sym4(i))
        recsym(i) = (-1+1i);
    end
    if(sym2(i)>=sym3(i) && sym2(i)>=sym4(i) && sym2(i)>=sym1(i))
        recsym(i) = (1+1i);
    end
    if(sym3(i)>=sym4(i) && sym3(i)>=sym1(i) && sym3(i)>=sym2(i))
        recsym(i) = (1-1i);
    end
    if(sym4(i)>=sym1(i) && sym4(i)>=sym2(i) && sym4(i)>=sym3(i))
        recsym(i) = (-1-1i);
    end
end


%SER

for j=1:(length(Tx))
    if(real(Tx(j)) < 0)
        if(imag(Tx(j)) < 0)
            Tx(j) = (-1 -1i);
        else
            Tx(j) = (-1 + 1i);
        end
    else
        if(imag(Tx(j)) > 0)
            Tx(j) = (1 + 1i);
        else
            Tx(j) = (1 - 1i);
        end
    end
end

count=0;
for i=1:length(Tx)
    if(Tx(i) ~= recsym(i))
        count = count+1;
    end
end
count




figure()
plot(real(recsym),imag(recsym),'*')
xlim([-1.5 1.5])
ylim([-1.5 1.5])

